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Abstract 

A new release of the Monte Carlo program HERWIG (version 6.3) is now available. The 
main new features are new (MRST) built-in parton distribution functions, more SM gauge 
boson production processes, 2^3 MSSM Higgs production processes, an option to handle 
negative event weights, and an interface to the CIRCE beamstrahlung program. 
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1 Introduction 

The last major public version (6.2) of HERWIG was reported in detail in [|]]. In this note we 
describe the main modifications and new features included in intermediate versions and the 
latest public version, 6.3. 

Please refer to [|lj and to the present paper if using version 6.3 of the program. 



1.1 Availability 

The new program can be obtained from the following web site: 



tittp : //hepwww . rl . ac . uk/theory/seymour/herwig/ 
This will temporarily be mirrored at CERN for the next few weeks: 



tittp : //home . cern . ch/seymour/herwig/ 



1.2 Increased common block size 

To take account of the increased energy and complexity of interactions at LHC and future 
colliders, the default value of the parameter NMXHEP, which sets the array sizes in the standard 
/HEPEVT/ common block, has been increased to 4000. 



1.3 ISAWIG update 

ISAWIG is an interface which allows SUSY spectra and decay tables generated by ISAJET f2[ 
to be read into HERWIG. A new version (1.104) of ISAWIG which works with ISAJET7.51 and 
includes pions for AMSB models has been released. (Note: we have been unable to generate 
SUGRA point 6 with this version of ISAJET.) 

The previous version of ISAWIG, 1.103, fixed a bug in chargino and sfermion mixing matrices. 
The problem was in the conversion between the convention of ISAJET and Haber and Kane. 

For fuller details see the ISAWIG web site: 



http : //www-thphys . physics . ox . ac . uk/users/PeterRichardson/HERWIG/isawig . html 



2 New Parton Distribution Functions 

The default parton distributions in HERWIG were very old and did not include fits to any of 
the HERA data. In the past users could link to PDFLIB in order to use more recent PDFs. 
However as this no longer seems to be maintained several new PDFs have been included in the 
new version. 



NSTRU 



Description 



6 

7 



Central as and gluon leading-order fit of [gj 

Higher gluon leading-order fit of || 

Average of central and higher gluon leading-order fits of 



It should be noted that we have only added leading-order fits because the evolution al- 
gorithms in HERWIG, in particular the backward-evolution algorithm for initial-state parton 
showering, are only leading-order and therefore inconsistencies could occur with next-to- leading- 
order distributions. 

The new default structure function set NSTRU=8 is the average of two of the published fits 
@, because this has been found || to be closer to the central value of more recent next-to- 
leading-order fits. The other fits can then be used to assess the effects of varying the high-x 
gluon. 

3 New Processes 

3.1 New gauge boson production processes 

Standard Model gauge boson production gives rise to important backgrounds for new physics 
searches at LHC. Two new classes of processes are included. 



3.1.1 IPROC=2800-2825: Gauge boson pair production 

The code already included in HERWIG for e + e~ — > WW/ZZ |J was adapted for hadron-hadron 
collisions, including photons and the photon/Z interference for the resonant diagrams. 



IPROC 


Process 


2800 


WW production in hadron-hadron collisions 


2810 


ZZ production in hadron-hadron collisions (including photon terms) 


2815 


ZZ production in hadron-hadron collisions (Z only) 


2820 


WZ production in hadron-hadron collisions (including photon terms) 


2825 


WZ production in hadron-hadron collisions (Z only) 



All of these processes use a cut EMMIN (default value 20 GeV) on the mass of the gauge 
bosons produced. The cut PTMIN (default 10 GeV) on the transverse momentum of the bosons 
is also used. Both these cuts should not be taken to zero simultaneously if photon terms are 
included. The phase space for these processes contains a number of peaks and it was therefore 
necessary to use an adaptive multi-channel phase-space integration method which is described 
below. 

A number of new subroutines, given in Table |T|, were added for these processes. 

3.1.2 IPROC=2900-2916: ZQQ production 



IPROC 


Process 


2900+IQ 


99 + QQ- 


-> ZQQ for massless Q and Q (IQ=1. . . 6 for Q = d . . .t) 


2910+IQ 


gg + qq- 


■* ZQQ, for massive Q and Q (IQ=1. . . 6 for Q = d . . .t) 



The matrix elements of || were used for the massless case and an independent calculation, 
using the approach of Q] , which was checked both for gauge invariance and against the massless 
case for the massive result. 

In both cases the decay of the Z is fully included and is selected using M0DB0S. PTMIN 
controls the minimum transverse momentum of the outgoing quarks. As with gauge boson pair 
production it was necessary to use an optimized multi-channel phase-space integrator which is 
described below. 

A number of new subroutines, given in Table [I], were added for these processes. 



3.1.3 Optimized phase space 

The phase space for both gauge boson pair production and ZQQ is complicated, as these 
processes are both treated as 2 — > 4 processes. In order to obtain a reasonable efficiency it was 
necessary to adopt a multi-channel approach based on that described in || [| . 

In each case a number of different channels are included which attempt to map the phase 
space for the different processes. The default weights for these different channels have been cho- 
sen to optimize the efficiency for the Tevatron and LHC; a choice of which to use is made based 
on the beam energy. However, the choice is affected by the phase space cuts applied. Therefore 
if these are significantly altered the weights for the different channels need re-optimizing. 

This is controlled by the new variable 0PTM (default .FALSE.). If 0PTM=.TRUE., before 
performing the initial search for the maximum weight HERWIG will attempt to optimize the 
efficiency using the procedure suggested in M. 



This is done by generating IOPSTP (default 5) iterations of IOPSH (default 100000) events. 
The choice of IOPSTP and IOPSH is a compromise between run time and accuracy. The value 
of IOPSH should not be significantly reduced because the procedure attempts to minimize the 
error on the Monte Carlo evaluation of the cross section and if IOPSH is small the error on the 
error can be significant. If you need to re-optimize the weights we would recommend a long 
run just to optimize the weights which can then be used in all the runs to generate events. The 
new subroutine HWIPHS was added to initialize the phase space. 



3.2 New 2 — ► 3 MSSM Higgs production processes 

A large number of final states involving the production of both neutral and charged Higgs 
bosons of the Minimal Supersymmetric Standard Model (MSSM) have been made available 
in version 6.3. They all proceed via 2 — >■ 3 body hard scattering subprocesses. They are 
listed below, with corresponding process numbers (IQ and ISQ are as detailed in the following 
subsections). Further details of their implementation can be found in [|Kj . 

The new subroutines introduced to administer the following processes are HWHIBQ and 
HWH2BH for IPR0C=3500, plus HWHISQ and HWH2SH for IPR0C=3100, 3200. In addition, HWHIGQ 
has been modified to accommodate the IPR0C=3800 series. Finally, the matrix element 
NME=200, describing the 1 — i- 3 body heavy-quark decays via a virtual H ± boson is now 
available and a new function, HWDHWT, has been introduced to this end. This can be used to 
emulate e.g. t —> 6iJ + (— > //') decays. 



IPRDC 


Process 


3100+ISQ 


99 /qq -* qq*H ± (ISQ=IPR0C-3100 as from first table below) 


3200+ISQ 


gg/qq -* qq*h, H, A (lSQ=IPR0C-3200 as from second table below) 


3500 


bq — ■> bq' H ± + ch. conj. 


3710 


qq -> q'q'h 


3720 


qq -> q'q'H 


3810+IQ 


gg + qq ~~ * QQh (all q flavours in s-channel, IQ as usual for Q flavour) 


3820+IQ 


gg + qq^QQHC) 


3830+IQ 


gg + qq^QQAC) 


3839 


gg + qq — > btH + + ch. conjg. (all q flavours in s-channel) 


3840+IQ 


gg ~^ QQh (IQ as above) 


3850+IQ 


gg^QQH (") 


3860+IQ 


gg^QQAC) 


3869 


gg — > btH + + ch. conjg. 


3870+IQ 


qq ~^ QQh (all q flavours in s-channel, IQ as above) 


3880+IQ 


qq^QQH (") 


3890+IQ 


qq^QQAC) 


3899 


qq — ► btH + + ch. conjg. (all q flavours in s-channel) 



3.2.1 IPROC=3110— 3178: charged Higgs plus squark pair production 

The production of charged Higgs bosons of the MSSM in association with squark pairs, of 
bottom and top flavours only, is implemented via the 3100 series of IPR0C numbers, as follows. 



IPROC 


partons - 


-> spartons 


Higgs 


3110 


gg + qq - 


-> ft?/ 


ff± 


3111 


99 + QQ ~ 


-» 6it* 


#+ 


3112 


99 + QQ ~ 


-> M2 


#+ 


3113 


gg + qq - 


-> Ml 


#+ 


3114 


gg + qq - 


■* M2 


H+ 


3115 


gg + qq - 


-» m 


H- 


3116 


gg + qq - 


-» £i&* 


H~ 


3117 


gg + qq - 


-> t 2 6* 


H~ 


3118 


gg + qq - 


"> «2&2 


H~ 



Add 30(60) to IPROC for gg(qq)-only initiated processes 



Their phenomenological relevance has been discussed in |[J 



3.2.2 IPROC=3210— 3298: neutral Higgs plus squark pair production 

The production of neutral Higgs bosons of the MSSM in association with squark pairs, of 
bottom and top flavours only, is implemented via the 3200 series of IPROC numbers, as follows. 



IPROC 


partons - 


-> spartons 


Higgs 


3210(3220) [3230] 


gg + qq - 


■* ¥j 


/i(iy)[A] 


3211(3221)[3231] 


gg + qq - 


-> Mi 


/i(#)[A] 


3212(3222) [3232] 


gg + qq - 


- hb* 2 


h{H)[A] 


3213(3223) [3233] 


gg + qq - 


- Mi 


h{H)[A] 


3214(3224) [3234] 


gg + qq - 


-> b 2 b* 


h{H)[A) 


3215(3225) [3235] 


gg + qq - 


-> tit* 


h(H)[A] 


3216(3226) [3236] 


gg + qq - 


-> ii<2 


h{H)[A] 


3217(3227) [3237] 


gg + qq - 


-> *2«l 


h{H)[A] 


3218(3228) [3238] 


gg + qq - 


■* ^2*2 


h(H)[A] 



Add 30(60) to IPROC for gg(qq)-on\y initiated processes 



Their phenomenological relevance has been discussed in [11, 12 



3.2.3 IPROC=3500: charged Higgs boson from 6g-initiated processes 



This process is relevant for charged Higgs scalar production at large tan/3 values, see [13 



3.2.4 IPROC=3710— 3720: neutral Higgs production via weak boson fusion 

These processes are the MSSM counterparts of the SM process of weak vector- vector fusion in 
hadronic collisions (IPR0C=1900). They are computed using the same subroutines and setting 
the $°VV couplings to sin(/3 - a) for $° = h° (IPR0C=3710) and to cos(/3 - a) for $° = H° 
(IPR0C=3720), respectively, where V = W , Z°. (There is no A°VV coupling at tree level.) 

These reactions are two of the major direct production channels of neutral CP-even Higgs 
bosons of the MSSM at hadron colliders, such as the LHC (see e.g. ||14|| ). 



3.2.5 IPROC=3811— 3899: Higgs boson plus heavy quark pair production 

For the case of neutral Higgs states, IPR0C=3810+IQ corresponds to h° production, 
IPR0C=3820+IQ to H° and IPR0C=3830+IQ to A . (For the last case, the variable PARITY 
is set to — 1.) Note also the production of charged Higgs states, via IPR0C=3839, 3869 and 
3899, in association with pairs of top-bottom quarks. 

In the usage of the IPROC numbers corresponding to neutral Higgs states, when 6-quarks 
are involved in gg-fusion modes (IPR0C(+30)=3845, 3855 or 3865), the user should take care 
to avoid double- counting the chosen process with the corresponding 2 — ► 1 and 2 — > 2 cases of 
IPR0C=3610-3630 and IPR0C=3410-3430 initiated by quark-antiquark annihilation, i.e. bb — > 
Higgs, and (anti)quark-gluon scattering, i.e. bg — > b Higgs, respectively: see |D|. Similar 
arguments hold for the charged Higgs states, as the gg- induced process (IPR0C(+30)=3869) is 
an alternative implementation of IPR0C=3450 fll6"| . 

The associate production of neutral Higgs bosons (both CP-even and CP-odd) of the MSSM 
with heavy QQ pairs (Q = b and t) is of extreme phenomenological relevance as a discovery 
mode of Higgs scalars, both at the Tevatron (Run 2) and the LHC (see e.g. [Ll, [IT]]), as is the 



case of the charged Higgs channel [16, 18 



4 Negative Weights Option 

A Monte Carlo program generates N w weights {wi} such that the estimated cross section is 

I N w 

(J =TrY. w ^ = w ■ (!) 



The corresponding error depends on the width of the weight distribution: 

Sa 1 Sw r 



■'mis 



o vAL w 



(2) 



If only positive weights are generated, and there exists a maximum weight u> max , then un- 
weighted events can be generated by 'hit and miss': w\ = or 1 with probability P{w' i = 1) = 
Wi/w max . Then 

8a /w ma .x - w lw mxx — w 1 , - 



r^ 



a V N W W V ^Wmax V^e 

where N e = N w w/w maiX is the number of events generated. The time needed (especially for 
detector simulation) depends mainly on the number of events. Hence the inefficiency of 'hit 
and miss' is not necessarily a disaster. This is the usual approach adopted in Monte Carlo 
event generators. 

Negative weights can be generated by subtraction procedures for matrix element corrections. 
These are not a problem of principle but prevent naive 'hit and miss'. To generalize 'hit and 
miss', one can generate unweighted events (w^ = 1) and 'antievents' (w't = —1) with 

sigri(w-) = sign(w 4 ) , (4) 

P(w'i = ±l) = H/M max . (5) 

Then 

5a_ _ MHmax- W 2 1 \w\ 

a ~ \ N w w 2 ~ y/W e w [ ' 



where N P = N„ 



\w\ \w\ 



is now the total number of events+ antievents generated. Again, the 



time needed is almost proportional to N e , so this is tolerable as long as \w\ 
section after any cuts that may be applied is 



w. 



Or 



>+-»- 



The cross 



(7) 



where N + events and N_ antievents pass the cuts. 

To allow for the possibility of negative weights, a new logical parameter NEGWTS has been 
introduced. The default (NEGWTS= . FALSE . ) is as before: negative weights are forbidden. If one 
is detected, a non-fatal warning is issued and the event weight is set to zero. 

If NEGWTS= . TRUE . , negative weights are allowed. Statistics are computed and printed ac- 
cordingly. If unweighted events are requested (N0WGT= . TRUE . ) , the initial search stores the 
maximum and mean absolute weights, \w | max and |iu|. Events and antievents are selected ac- 
cording to Eqs. (fUD and EVWGT is reset to |«;|sign(iUj), so that the numerator in Eq. (|7|) is the 
sum of EVWGTs for contributing (anti)events. 



5 CIRCE Interface 



Simulation of beamstrahlung is now available, via an interface to the CIRCE program |Lj|. It 
is implemented as a modification of the standard bremsstrahlung implementation, so is only 
available for processes for which that is available. The produced radiation is treated in the 
collinear limit, i.e. COLISR and the option WWA in routine HWEGAM are forcibly set .TRUE. . 

For both types of distribution function, f e / e and / 7 / e > we use a simplifying approximation: 
where the correct result should use the convolution of beamstrahlung and bremsstrahlung, we 
actually use the sum of the two. For the photon distribution this is quite reasonable, but for 
the electron it is perhaps more questionable. It boils down to the replacement: 

Je/e\X) / Je/e,breml Z J Je/e,beam I ) * Je/e,brem\X J T Je/e,bea,m\X) "v^ X). [Of 

Jx Z ^ ' ^ Z' 

This is good to the extent that both distributions are dominated by the x — > 1 region. A small 
utility program can be obtained from the HERWIG web page to test this approximation. For 
example, for TESLA running at 500 GeV, the integral over f e / e (x) for large x, relevant for 
threshold shapes, is accurate to 1% for x > 0.95 and to 10% for x > 0.99, and the value of 
/e/eUfJoo) 2 ), re l evan t f° r radiative return to the Z, is accurate to 2%. 

In order for the distribution to be probabilistic, the coefficient of the delta function has 
to remain positive. If it is not, HERWIG will terminate. This can be prevented by adjusting 
the ZMXISR parameter: the default value, ZMXISR=0. 999999 is too large and HERWIG will fail, 
but reducing it to ZMXISR=0. 99999 works fine for all the parameter sets currently available. It 
is straightforward to check that ZMXISR values in this range do not have an influence on any 
physical observables, by rerunning with it further reduced. 

The interface is controlled by five new variables, which are given in the following table 
together with their default values: 



CIRC0P 





CIRC AC 


2 


CIRCVR 


7 


CIRCRV 


9999 12 31 


CIRCCH 






CIRCOP is the main control option: CIRC0P=0 switches off beamstrahlung and uses stan- 
dard HERWIG; CIRC0P=1 switches to collinear kinematics, but still leaves beamstrahlung 
switched off; CIRC0P=2 uses only beamstrahlung; and CIRC0P=3 uses both beamstrahlung 
and bremsstrahlung. CIRC0P=0 and CIRC0P=3 should therefore be regarded as 'off' and 'on', 
respectively, with the other two options mainly for cross-checking purposes. The variables 
CIRCAC, CIRCVR, CIRCRV and CIRCCH are simply passed to CIRCE as its input variables ace, 
ver, rev and chat, as described in its documentation. The default values correspond to the 
most up-to-date revision of version 7 of the TESLA parametrization, with minimal output. 
CIRCE can be obtained from 



tittp: //heplix. ikp.physik.tu-darmstadt .de/nlc/beam.html 



6 Minor Changes and Bug Fixes 

• In the INCLUDE file, common blocks have been regrouped so that all commons that were 
present in version 6.1 are unchanged. All new variables since version 6.1 are in /HWGRAV/, 
/HWPMRS/, /HWCIRC/, /HW6202/, /HW6203/ or /HW6300/. From now on all common blocks 
will be frozen and new variables introduced in version x.yyy will be put in a new common 
block /HWxyyy/. 

• In HWHDYP lepton mass effects are corrected. 

• In HWHIGS a factor of 2 excess in the cross section for Higgs production by gluon or quark 
fusion is corrected. 

• In HWUDKL safety against numerical overflows for long-lived particles is improved. 

• In HWUMBW Breit-Wigner smearing of W and Z boson masses is enabled. 

• A work-around to overcome problems with the new DEC UNIX Fortran compiler has 
been implemented. 

• In IPR0C=3610-30, top quark decays now take place. Previously tops produced in these 
processes were missed by the decayer and formed top hadrons. In addition, changes in 
the value of IPROC during generation of these processes have been eliminated. 

• Gluon spin correlations have been corrected in gg — > A (IPR0C=3630) and implemented 
in gg -> Graviton (IPR0C=4200 etc.). 

• Underlying event suppression by IPROC^IPROC+10000 is now enabled in all SUSY pro- 
cesses. 

• CKM matrix element bugs in W+Higgs and heavy quark+Higgs production 
(IPROC/ 100=33,34) have been corrected. 

• The default maximum number of errors in event generation is now set to 
MAX(10,MAXEV/100). 

• A dummy time subroutine TIMEL has been included. It should be deleted and replaced 
by a system or CERN Library routine giving TRES= CPU time remaining (seconds) if 
this is needed, e.g. to terminate batch jobs. The dummy returns TRES=10 10 . 



Routine 



Description 



CIRCEE Dummy CIRCE routine - delete if using CIRCE 

CIRCES Dummy CIRCE routine - delete if using CIRCE 

CIRCGG Dummy CIRCE routine - delete if using CIRCE 

HWDHWT Subroutine for top decay via a virtual H ± 

HWH2BH Matrix element for H^ production via 6g-fusion 

HWH2DD Function to return the D function of || 

HWH2F1 Subroutine to return the F function of |7| for a fixed first momentum 

HWH2F2 Subroutine to return the F function of |7j for a fixed second momentum 

HWH2F3 Subroutine to return the F function of |7| for all first and second momenta 

HWH2M0 Subroutine to compute the massless matrix element for QQZ 

HWH2MQ Subroutine to compute the massive matrix element for QQZ 

HWH2PS Subroutine to perform the phase-space for Z+two jets 

HWH2P1 Subroutine to select quark masses for HWH2PS 

HWH2P2 Subroutine to select quark masses for HWH2PS 

HWH2SH Matrix element for squark pair plus Higgs production 

HWH2SS Subroutine to return the S function of [pj 

HWH2T1 Function to return the T\ function of || 

HWH2T2 Function to return the T 2 function of [EJ 

HWH2T3 Function to return the T 3 function of [BJ 

HWH2T4 Function to return the T4 function of || 

HWH2T5 Function to return the T5 function of || 

HWH2T6 Function to return the T 6 function of || 

HWH2T7 Function to return the T 7 function of [EJ 

HWH2T8 Function to return the T 8 function of [BJ 

HWH2T9 Function to return the Tg function of [EJ 

HWH2T0 Function to return the Tio function of 

HWHDYQ Subroutine for QQZ 

HWHGBP Main routine for gauge boson pair production in hadron-hadron 

HWHGBS Phase space for gauge boson pair production in hadron-hadron 

HWHGB1 Selects gauge boson mass for HWHGBS 

HWHGB2 Matrix element for WW in hadron-hadron 

HWHGB3 Matrix element for ZZ in hadron-hadron 

HWHGB4 Matrix element for WZ in hadron-hadron 

HWHGB5 Selects t and u for HWHGBS 

HWHIBQ Subroutine for H^ production via 6g-fusion 

HWHISQ Subroutine for squark pair plus Higgs production 

HWHV2J Master subroutine for all gauge boson + two jet processes 

HWIPHS Subroutine to initialize the optimized phase space 

HWSMRS Subroutine for MRST PDFs 

TIMEL Dummy time subroutine - see Sect. 6 



Table 1: New subroutines for version 6.3. 
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